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High-accuracy numerical simulation of black-hole binaries: Computation of the 
gravitational-wave energy flux and comparisons with post-Newtonian approximants 
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Expressions for the gravitational wave (GW) energy flux and center-of-mass energy of a compact 
binary are integral building blocks of post-Newtonian (PN) waveforms. In this paper, we compute 
the GW energy flux and GW frequency derivative from a highly accurate numerical simulation of 
an equal-mass, non-spinning black hole binary. We also estimate the (derivative of the) center- 
of-mass energy from the simulation by assuming energy balance. We compare these quantities 
with the predictions of various PN approximants (adiabatic Taylor and Pade models; non-adiabatic 
effective-one-body (EOB) models). We find that Pade summation of the energy flux does not 
accelerate the convergence of the flux series; nevertheless, the Pade flux is markedly closer to the 
numerical result for the whole range of the simulation (about 30 GW cycles). Taylor and Pade 
models overestimate the increase in flux and frequency derivative close to merger, whereas EOB 
models reproduce more faithfully the shape of and are closer to the numerical flux, frequency 
derivative and derivative of energy. We also compare the GW phase of the numerical simulation 
with Pade and EOB models. Matching numerical and untuned 3.5 PN order waveforms, we find that 
the phase difference accumulated until Mui = 0.1 is -0.12 radians for Pade approximants, and 0.50 
(0.45) radians for an EOB approximant with Keplerian (non-Keplerian) flux. We fit free parameters 
within the EOB models to minimize the phase difference, and confirm the presence of degeneracies 
among these parameters. By tuning the pseudo 4PN order coefficients in the radial potential or 
in the flux, or, if present, the location of the pole in the flux, we find that the accumulated phase 
difference at Mu = 0.1 can be reduced — if desired — to much less than the estimated numerical 
phase error (0.02 radians). 

PACS numbers: 04.25.D-, 04.25.dg, 04.25.Nx, 04.30.-w 



I. INTRODUCTION 

The first-generation interferometric gravitational wave 
(GW) detectors, such as LIGO [1, 2], GEO600 [3] and 
Virgo [4, 5], are now operating at or near their design 
sensitivities. One of the most promising sources for these 
detectors is the inspiral and merger of binary black holes 
(BBHs) with masses m 1 ~ m 2 ~ 10-20 M Q [6, 7]. A de- 
tailed and accurate understanding of the gravitational 
waves radiated as the black holes spiral towards each 
other will be crucial not only for the initial detection 
of such sources, but also for maximizing the information 
that can be obtained from signals once they are observed. 
Both the detection and subsequent analysis of gravita- 
tional waves from compact binaries depends crucially on 
our ability to build an accurate bank of templates, where 
each template is a theoretical model that accurately rep- 
resents the gravitational waveform from a binary that has 
a certain set of parameters (e.g., masses and spins). For 
detection, the technique of matched filtering is applied to 
noisy data to extract any signals that match members of 
the template bank. For analysis, the best-fit parameters 
are determined, most likely by an iterative process that 
involves constructing further templates to zero in on the 
best fit. 

When the black holes are far apart and moving slowly, 



the gravitational waveform (i.e., the template) can be 
accurately computed using a post-Newtonian (PN) ex- 
pansion. As the holes approach each other and their 
velocities increase, the post-Newtonian expansion is ex- 
pected to become less and less reliable. However, until 
recently there has been no independent way to deter- 
mine how close comparable-mass holes must be before 
PN methods become inaccurate. This has changed with 
recent advances in numerical relativity (NR), which make 
it possible for the first time to quantify the disagreement 
between PN predictions [8] and the true waveform [9-14] . 
In a previous paper [12], some of us described numerical 
simulations of 15 orbits of an equal-mass non-spinning 
binary black hole system. Gravitational waveforms from 
these simulations covering more than 30 GW cycles and 
ending about 1.5 GW cycles before merger, were com- 
pared with those from quasi-circular PN formulas for 
several time-domain Taylor approximants computed in 
the so-called adiabatic approximation. We found that 
there was excellent agreement (within 0.05 radians) in 
the GW phase between the numerical results and the PN 
waveforms over the first ~ 15 cycles, thus validating the 
numerical simulation and establishing a regime where PN 
theory is accurate. In the last 15 cycles to merger, how- 
ever, generic time-domain Taylor approximants build up 
phase differences of several radians. But, apparently by 
coincidence, one specific PN approximant, TaylorT4 at 



3.5PN order, agreed much better with the numerical sim- 
ulations, with accumulated phase differences of less than 
0.05 radians over the 30-cycle waveform. Simulations by 
Hannam et al. [14] for equal-mass, non-precessing spin- 
ning binaries confirm that this agreement in the non- 
spinning case is a coincidence: they find the phase dis- 
agreement between TaylorT4 and the numerical wave- 
form can be a radian or more as the spins of the black 
holes arc increased. 

To build a template bank to be used by ground-based 
GW detectors, one possibility would be to run a sep- 
arate numerical simulation for each template. This is 
not currently possible, however, due to the large com- 
putational cost per numerical waveform (on the order 
of a week for a single waveform) and the large number 
of templates needed to cover the parameter space, es- 
pecially when spins are present. A more realistic possi- 
bility is to perform a small number of simulations and 
develop an analytic template family (i.e., a fitting for- 
mula) which interpolates the parameter space between 
the simulations [15-20]. 

Before the NR breakthrough several analytic prescrip- 
tions were proposed to address the loss of accuracy 
of the adiabatic Taylor approximants. Damour, Iyer 
and Sathyaprakash [21] introduced the Pade summation 
of the PN center-of-mass energy and gravitational en- 
ergy flux in order to produce a series of Pade approx- 
imants for the waveforms in the adiabatic. Buonanno 
and Damour [22-25] introduced the effective-one-body 
(EOB) approach which gives an analytic description of 
the motion and radiation beyond the adiabatic approx- 
imation of the binary system through inspiral, merger, 
and ringdown. The EOB approach also employs the 
Pade summation of the energy flux and of some cru- 
cial ingredients, such as the radial potential entering 
the conservative dynamics. So far, the EOB waveforms 
have been compared with several numerical waveforms 
of non-spinning binary black holes [9, 15, 16, 18-20]. 
Buonanno et al. [16] showed that by using three quasi- 
normal modes [9] and by tuning the pseudo 4PN order 
coefficient [26] in the EOB radial potential to a specific 
value, the phase difference accumulated by the end of the 
ringdown phase can be reduced to ~ 0.19-0.50 radians, 
depending on the mass ratio and the number of multi- 
pole moments included in the waveform. Those results 
were obtained using waveforms with 5-16 GW cycles and 
mass ratios 1 : 4, 1 : 2, 2 : 3 and 1 : 1. In Refs. [18-20] 
the authors introduced other improvements in the EOB 
approach, in part obtained by tuning the test-mass limit 
results [27] - - for example Pade summation of the PN 
amplitude corrections in the inspiral waveform; ringdown 
matching over an interval instead of a point; inclusion 
of non-circular terms in the tangential damping force; 
use of five quasi-normal modes. They found that the 
phase differences accumulated by the end of the inspiral 
(ringdown) can be reduced to ±0.001 (±0.03) radians for 
equal-mass binaries [18, 19] and to ±0.05 radians for bi- 
naries with mass ratio 1 : 2 [20]. Note that these phase 



differences arc smaller than the numerical errors in the 
simulations. 

The energy flux and the center-of-mass energy are two 
fundamental quantities of the binary dynamics and cru- 
cial ingredients in building GW templates. In this pa- 
per we extract these quantities, and compare the results 
from our numerical inspiral simulation [12] with PN re- 
sults in both their Taylor-expanded and summed (Pade 
and EOB) forms. The agreement between the numerical 
and analytical results for the energy flux and the center- 
of-mass energy is a further validation of the numerical 
simulation. It also allows us to study whether or not the 
agreement of the phase evolution of PN and numerical 
waveforms is accidental. In addition, we compute wave- 
forms based on adiabatic Pade and non-adiabatic EOB 
approximants in their untuned form (i.e., without in- 
troducing fitting coefficients) and study their agreement 
with our numerical simulations. 

We try to understand whether these approximants can 
reproduce features of the numerical simulations that can 
be exploited to develop a faithful analytic template fam- 
ily. By introducing unknown higher-order PN coefficients 
into the dynamics and tuning them to the numerical data, 
we investigate how to improve the agreement with the nu- 
merical results. Although our study only examines non- 
spinning, equal-mass binary black holes, by combining it 
with other studies [15-20] one can already pinpoint which 
parameters are degenerate and which have the largest ef- 
fect on the waveforms. This is particularly relevant dur- 
ing the last stages of inspiral and plunge. The overall 
methodology can be extended to a larger region of the 
parameter space. We will defer to a future paper a com- 
plete study of the flexibility of the EOB approach with 
the extension of our numerical waveform through merger 
and ringdown. 

This paper is organized as follows: Section II gives 
a quick review of the numerical simulations presented 
in [12], and then presents the computation of the GW 
energy flux from the simulation. In Sec. Ill we summa- 
rize the PN approximants that will be compared to the 
numerical simulation. In Sec. IV, wc compare the GW 
energy flux for the various PN approximants with numer- 
ical results and explore the possibility of improving the 
agreement with the numerical flux by adding phenomeno- 
logical parameters [15, 16, 18-20]. In Sec. V, we examine 
the evolution of the center-of-mass energy for the various 
PN approximants and compare to the numerical results 
assuming balance between the change in the center-of- 
mass energy and the energy carried from the system by 
the gravitational waves. In Sec. VI we compare wave- 
forms constructed from the Pade and EOB approximants 
with our numerical results, and study how to improve the 
agreement by exploiting the flexibility of the EOB model 
(i.e., by fitting free parameters of the EOB model). Fi- 
nally, we present some concluding remarks in Sec. VII. 
In the Appendix we review the performance of the Pade 
summation of the Taylor series of the energy flux in the 
test particle limit. 
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FIG. 1: Some aspects of the numerical simulation. From 
top panel to bottom: the leading mode /122; the two next 
largest modes, /144 and /132 (smallest); the frequency of /122 
[see Eq. (5)]. 



II. COMPUTATION OF THE NUMERICAL 
GRAVITATIONAL- WAVE ENERGY FLUX 

A. Overview and Definitions 

The data used in this paper is the same as that de- 
scribed in Sec. II of Boyle et al. [12]. The simulation is 
a 16-orbit inspiral, with very low spin and eccentricity. 
Figure 1 presents a view of some relevant quantities of 
that simulation. 

The Newman- Penrose scalar ^4, defined using a 
coordinate-based tetrad, is extracted from the simulation 
at several extraction radii and expanded in spin- weighted 
spherical harmonics, 
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Then \I/ 4 m (£, r) is extrapolated to infinite extraction ra- 
dius using an n-th order polynomial in 1/r, where typ- 
ically n = 3. This results in the asymptotic field 
r\& 4 m (t — r* ) as function of retarded time 1 t — r* . 

Gravitational radiation may also be expressed via 
the standard metric-perturbation quantities h+ and h x , 
which we similarly write in terms of spin- weighted spher- 



1 See Sec. II F of Ref. [12] for a precise definition of r* and a 
description of the extrapolation. 
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For linear perturbations around Minkowski space, 
* 4 m (i — r*) — hi m (t — r*). In particular, this relation 
should be true for the waveforms we have extrapolated 
to infinity. 

However, to compute the energy flux we do not need 
to determine h; we need only its time derivative h. The 
energy flux depends on the spin-weighted spherical har- 
monic coefficients of the time derivative h via 
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as discussed 



We obtain hi m by time-integration of ^4 
in detail below. 

Finally, we define gravitational wave phase and fre- 
quency in two ways — one based on ^4 2 , and one based 
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In both cases, we define the arg function to be the usual 
function, with discontinuities of 2-77 removed. Many PN 
formulae (see Sec. Ill) involve yet another frequency and 
phase: the orbital phase 4> and orbital frequency f2. Al- 
though the three frequencies satisfy u> « w « 2(7, the 
slight differences between different frequencies are sig- 
nificant at the level of precision of our comparison (see 
Fig. 6 below), so it is important to distinguish carefully 
between them. 

When discussing our numerical solution, we write all 
dimensionful quantities in terms of the mass scale M, 
which we choose to be the sum of the irreducible masses 
of the two black holes. 2 



B. Calculation of h 

The energy flux depends on the spin- weighted spherical 
harmonic coefficients of h via Eq. (3). We therefore need 
to perform one time integration on ^ 4 m : 



hi m (t) 
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This integration is performed for each mode (l,m) sep- 
arately and requires the choice of two integration con- 
stants, which arc contained in the complex number Hi m . 



This quantity was denoted by m in Ref. [12] 



Ideally. Hi m should be chosen such that hi m —> for 
t — > — oo. Because our numerical simulations do not ex- 
tend into the distant past, this prescription cannot be 
implemented. Rather, we make use of the approximation 
that the real and imaginary parts of hi m should oscillate 
symmetrically around zero. 

Let us consider a pure sine/cosine wave, with constant 
amplitude and phase: 
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A[cos(u>t) + i sin(wt)] , (7) 

A 

- [sin(wi) - i cos(wi)] + ff cx , (8) 



where the superscript 'ex' stands for example. The am- 
plitude is given by 

A 2 A 
1^x12 = + 2_[R ej ff ex sin(^)-Imiy ex cos(wi)] + |H ex | 2 . 

LU 2 UJ 

(9) 
Only for the correct choice of integration constants, 
H ex = 0, is the amplitude \h ox \ constant. 

Therefore, we propose to determine the integration 
constants Hi m in Eq. (6) by minimizing the time deriva- 
tive of the amplitude over the entire waveform. In par- 
ticular we minimize 
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From this minimization principle it follows that Hi m is 
determined by the linear system 

ReH (Re^^fdt + ImH /Re^Im^d* 



(Re^^Re/io + Re* 4 Im<I' 4 Im/2r 



RcH Rc^ 4 Im^ 4 dt + ImH (Imf 4 ydt 



(Im* 4 ) 2 Im/i + RevI' 4 Im* 4 Rc/io 



dt. 
(11a) 

dt. 

(lib) 



Here, we have suppressed the indices Im for clarity, 
all integrals are definite integrals from t\ to t 2 , and 
ho(t) = J ^> 4 (t')dt'. For a given integration interval 
[ii,^]) Eqs. (11) provide a deterministic procedure to 
determine the integration constants Hi m . We note that 
there have been several earlier proposals to fix integra- 
tion constants [20, 28-31]. While we have not tested 
those proposals, we point out that Eqs. (11) allow for 
very accurate determination of the integration constants 
and one can easily obtain an error estimate, as we discuss 
in the next subsection. 



C. Uncertainties in numerical quantities 

Because the amplitude and frequency of the wave- 
form are not constant, this procedure is imperfect, and 



the result depends somewhat on the chosen values of 
t\ and t 2 . To estimate the residual uncertainty in H 
due to this choice, we select nine different values for t\ 
and eleven values for t 2 : h = 200M, 220M, . . . , 360M; 
t 2 = 2000M, 2100M, . . . , 3000M. The values of h vary 
over roughly one GW cycle and test the sensitivity to 
the GW phase at the start of the integration interval; 
the values of t 2 are designed to test the dependence 
on the amplitude at the end of the integration interval. 
For t 2 > 3000M we find that the errors in our proce- 
dure rapidly increase for several reasons: (a) the min- 
imization principle is based on the approximation that 
the amplitude is constant; this approximation becomes 
worse toward merger; (b) Zj TO in Eq. (10) weights abso- 
lute changes in \h\, not relative ones; close to merger, 
the amplitude becomes so large that it dominates X; m ; 
and (c) the integration constants shift the waveform hi m 
vertically, and we are trying to determine the particular 
vertical shift such that hi m is centered around zero. De- 
termination of such an offset is most accurate in a regime 
where the oscillations are small, i.e., at early times. 

For each of these 99 integration intervals, we compute 
integration constants using Eqs. (11) for the three dom- 
inant modes, h 22 , h± 4 and /132, and we compute F(t) 
from Eq. (3) using only these modes and we compute 
vj(t). (We will show below that the contributions of 
other modes are far below our numerical errors on the 
flux.) We average the 99 functions F(t) and w(t) and 
then use a parametric plot of F(t) versus w{t) in our 
comparisons presented below. The variation in these 99 
values yields an uncertainty in F due to the choice of 
integration constants. 

The lower panel of Fig. 2 shows the variation in flux 
from the 99 different integration intervals. We find that 
the maximum deviation can be well approximated by 
m&x\5F\/F = 1.5 x 10" 5 (Mg7)- 3 / 2 (see the solid line 
in lower panel of Fig. 2). The average F computed from 
all 99 intervals [£1,^2] will have a smaller error. Inspec- 
tion of the lower panel of Fig. 2 reveals that the 6F/F 
curves fall into 11 groups, corresponding to the 11 val- 
ues of t 2 . Assuming that SF between these groups is 
randomly distributed, the error of the average will be re- 
duced by a factor vTT, i.e., SF/F = 5 x 10~ 6 (.MV)- 3 / 2 . 
This error is indicated as the grey shaded area in the 
upper panel of Fig. 2. 

The upper panel of Fig. 2 plots the relative change in 
F(w) for several changes in our numerical simulation: 
(a) Computing the flux from a run with lower resolu- 
tion (0030c/N5 in the language of Boyle et al. [12]); (b) 
using a different set of extraction radii for the extrac- 
tion of the gravitational wave; (c) increasing the polyno- 
mial order of extrapolation of ^ 4 to infinite extraction 
radius from n = 3 to n = 4; and (d) computing the flux 
from a separate evolution with a different outer bound- 
ary radius (0030c-2/N6). At low frequencies, the error 
is dominated by extrapolation to infinite radius and is 
a few tenths of a percent; at intermediate frequencies, 
0.055 < Mw < 0.083, all errors are smaller than 0.1 per- 
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FIG. 2: Lower panel: Relative difference between flux F(w) 
computed with 99 different intervals [ti,*2] and the average 
of these. Upper panel: Relative change in the flux F(vu) 
under various changes to the numerical simulation. The grey 
area in the upper panel indicates the uncertainty due to the 
choice of integration constants, which is always dominated by 
numerical error. The dashed line in the upper panel is our 
final error estimate, which we plot in later figures. 
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FIG. 3: Contributions of various (Z,m)-modes to the total 
numerical gravitational wave flux. Upper panel: plotted 
as a function of time. Lower panel: Plotted as a function 
of frequency Mm, The lower panel also contains the error 
estimate derived in Fig. 2. 
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FIG. 4: Lower panel: Difference between frequency deriva- 
tive vu computed with 99 different intervals [ti , ta] and the 
average of these. Upper panel: Change in the frequency 
derivative zb under various changes to the numerical simula- 
tion. The grey area in the upper panel indicates the uncer- 
tainty due to choice of integration constants, which dominates 
the overall uncertainty for low frequencies. The dashed line 
in the upper panel is our final error estimate, which we plot 
in later figures. 



cent. At frequency Mm s=s 0.084 we change the gauge 
conditions in the evolutions to allow wave-escorting; this 
introduces high-frequency features, which are small when 
extrapolation order n = 3 is used, but which dominate 
for n = 4 extrapolation. The numerical data we use in 
the PN comparisons below is extrapolated with n = 3, for 
which the features due to change of gauge are small, but 
nevertheless we will use conservative error bars encom- 
passing the n = 4 extrapolation as indicated in Fig. 2, i.e. 
a relative error of 0.2 per cent for Mm > 0.083. We find 
that the uncertainty in the flux due to numerical error 
in determining ^4 is always larger than the uncertainty 
due to the choice of integration constants. 

The contributions of the various (I, m)-modes to the 
total flux [see Eq. (3)] are plotted in Fig. 3. The top panel 
plots the flux as a function of time; the lower panel as a 
function of frequency Mm. The dashed line in the lower 
panel corresponds to the error estimate of Fig. 2. Because 
the modes (5, 4), (6, 6), and (8, 8) are significantly smaller 
than our error estimate, we do not include them in the 
present analysis. 

To estimate the uncertainty in m, we proceed in a sim- 
ilar fashion. Each one of the 99 different integration in- 
tervals yields an h^i from which we determine m. We 
average these to obtain the final m to be used in the post- 
Newtonian comparisons. The lower panel of Fig. 4 shows 
the variation in m between the 99 different integration 



intervals. We find that the maximum deviation can be 
well approximated by max |M 2 <5zi7| = 5 x 10~ 6 (Mn7)~ - 3 
(see the solid line in lower panel of Fig. 4). The average w 
computed from all 99 intervals [tijfo] will have a smaller 
error. Inspection of the lower panel of Fig. 4 reveals that 
the 8zu curves fall into 11 groups, corresponding to the 11 
values of £2- As for the case of 5F, if we assume that 8w 
between these groups is randomly distributed, then the 
error of the average will be reduced by a factor vll, i.e., 
M 2 8w = 1.5 x Kr 6 (Mn7)- - 3 . This error is indicated as 
the grey shaded area in the upper panel of Fig. 4. 

The upper panel of Fig. 4 plots also the change in w{w) 
for the same changes in our numerical simulation already 
discussed above. We find that at Mw < 0.083, the un- 
certainty in w is dominated by the choice of integration 
constants, whereas at higher frequencies the uncertainty 
is dominated by the numerical errors in the calculation 
of ^4. As discussed above, at frequency Mw « 0.084 we 
change the gauge conditions in the evolutions to allow 
wave-escorting; this introduces high-frequency features 
leading to more conservative error estimates. 

Note that w is a very steep function of w. While the 
absolute errors in 137 are roughly constant for our sim- 
ulation, the relative errors change significantly: 8w/w 
drops from about 10 per cent early in the run to about 
0.2 percent at late times. 

We also point out that the first 1000M of our simula- 
tion are contaminated by noise due to a pulse of "junk- 
radiation" at the start of the simulation. While this con- 
tamination is not apparent on a plot of the waveform as 
in Fig. 1, it nevertheless limits accurate PN-NR compar- 
isons to the region, t - r* > 1000M, i.e., Mw > 0.037. 



III. POST-NEWTONIAN APPROXIMANTS 

In this paper we will compare the numerical simulation 
to various approximants based on the PN expansion. The 
PN expansion is a slow-motion, weak-field approximation 
to general relativity with an expansion parameter e ~ 
(v/c) 2 ~ (GM/rc 2 ). For a binary system of two point 
masses mi and m^, v is the magnitude of the relative 
velocity, M is the total mass, and r is the separation. 
For a review of the PN expansion applied to gravitational 
radiation from inspir ailing compact binaries, see Ref. [8]. 

In Table I we summarize the PN-approximants that 
we use, and our notation. We shall use the PN ap- 
proximants in the so-called adiabatic approximation, 
both in the standard Taylor-expanded form (reviewed in 
Sec. Ill A) and in a form based on Pade summation (re- 
viewed in Sec. IIIB). In addition we shall use the non- 
adiabatic EOB model (reviewed in Sec. Ill C) in its orig- 
inal form [22-24] , as well as several variations that differ 
in the form of the radiation-reaction force [32-34] . After 
summarizing the various PN approximants in Sees. Ill A, 
IIIB, and IIIC, we describe how we construct the wave- 
form for these approximants in Sec. HID. 

In the adiabatic approximation the inspiral is modeled 



approximant 


notation 


see Eqs. 


adiabatic 


Keplerian 


Taylor (T-) 


F n /E p 


(19)/(14) 


yes 


yes 


Pade (P-) 


F™/El 


(39)/(33) 


yes 


yes 


EOB (E-) 


F™/H p 


(64)/(44) 


no 


yes 


EOB (E-) 


nK F™/H p 


(65)/(44) 


no 


no 


EOB (E-) 


F n /H p 


(69)/(44) 


no 


yes 


EOB (E-) 


nK F n /H p 


(70)/(44) 


no 


no 



TABLE I: Summary of PN-approximants. The T- 

approximants are always Taylor T4 [12] except in Fig. 16. 
The P-approximant in the second row was introduced in 
Refs. [21, 24, 32] and the original E-approximant in third 
row was introduced in Refs. [22-24]. The last three rows refer 
to three possible variations of E-approximants introduced in 
Refs. [32, 33]. In a few tests aimed at improving the closeness 
between numerical data and E-approximants, we vary « po i 
and treat the logarithms as constants when Pade summation 
to the flux is applied [18]. We shall denote this flux by F n . 
Finally, when using tuned PN-approximants with pseudo 4PN 
order terms in the flux, energy, or Hamiltonian, we denote the 
latter as pF, pE and pH. Note that if known test-mass limit 
coefficients in the flux are used, the latter is still denoted as 
F even at PN orders larger than 3.5PN. Finally, the values of 
Unoit, and v\ so used in the P-approximants F™ and nK F™ are 
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= 0.6907 and^C = 0.4456. 



as a quasi-stationary sequence of circular orbits. The 
evolution of the inspiral (and in particular of the orbital 
phase 4>) is completely determined by the energy-balance 
equation [8] 



dE{v n ) 
dt 



-F(y n ) ■ 



(12) 



This equation relates the time derivative of the center- 
of-mass energy E(vq) (which is conserved in absence of 
radiation reaction) to the gravitational wave energy flux 
F(vq). Both functions are known for quasicircular orbits 
as a PN expansion in the invariantly defined velocity 



v Q = (Mft) 1/3 



(13) 



where tt = $ is the orbital frequency (we use units such 
that G = c = l). 3 We will denote the Taylor-expanded 
flux (energy) by Fk (Ek) where k denotes the maximum 
power of va retained in the series. (Recall that k = 2N 
for an Ath order PN expansion.) We will denote the 
Padc-expandcd flux (energy) by F™ (E™) where m+n = 
k, with m and n denoting the order of the polynomial in 
the numerator and denominator, respectively. 



A. Adiabatic Taylor approximants 

For generic values of the symmetric mass ratio v = 
mim2 /M 2 , the center-of-mass energy is known through 



In Ref. [12] we used x ■■ 



as the expansion parameter. 



3PN order [35-39]. For circular orbits the Taylor PN- 
approximants (henceforth, T-approximants) to the en- 
ergy are given by 



Mv 



E 2 k{vn) = ^-«n ^£2i{y)' 



'Q i 



where the known coefficients are 
£o(v) = 1, 
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(18) 



The GW energy flux for arbitrary masses has been 
computed through 3.5PN order [40, 41]: 



F k (v n ) 
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where je is Euler's constant. Notice that starting at 3PN 
order (k = 6) logarithms enter the flux. 



and by applying Pade summation [42] to build succes- 
sive approximants to these two functions (henceforth P- 
approximants). Their motivation for introducing these 
new functions and using their P-approximants came from 
an examination of the behavior of the standard PN- 
cxpansion and the new P-approximants in the test-mass 
limit in which the exact gravitational energy flux is 
known numerically [43] , the PN expansion of the flux is 
known through 5.5PN order [44], and the center-of-mass 
energy is known analytically as 



E(v n ;v = 0) 



1 



2^ 



yr 



"3^ 



1 



(28) 



where /i = Mv is the reduced mass. 

DIS first observed that in the quantum two-body prob- 
lem the symmetric quantity 



El 



2mi m 2 



(29) 



(where the total relativistic energy E to t = E + M) , is the 
best energy function when treating the two-body prob- 
lem as an effective one-body problem in an external field. 
Because in the test-mass limit 



1 - 1v% 

e{vn;v = {)) = 



v/T 



3^ 



DIS defined the new energy function as 



e(va) 



1 



(30) 



(31) 



as this function has a simple pole singularity on the real 
axis in the test-mass limit, and DIS conjectured that such 
a pole would continue to exist in the comparable mass 
case. 4 The energy function E(vq) entering the balance 
equation (12) can be expressed in terms of e(i>n) as 



E(va) = {m 2 



2vM z 



s/l + e(v n ) - l] } 



1/2 



-M. 



(32) 

by combining Eqs. (29) and (31). [Note that the map be- 
tween the adiabatic functions e and E given by Eq. (32) 
is the same map found in the EOB model between the cf- 
fective Hamiltonian H cS and the real Hamiltonian H^\ 
as given by Eq. (44).] 

Finally, DIS proposed as approximants to the en- 
ergy function e(wn) the diagonal or subdiagonal P- 
approximants, depending on whether the PN order is 
even or odd. 5 Investigating the behavior of the P- 



B. Adiabatic Pade approximants 

1. Center-of-mass energy 

Damour, Iyer and Sathyaprakash [21] (henceforth DIS) 
proposed a new class of approximate waveforms con- 
structed by introducing new energy and flux functions 



4 A motivation for having using Eq. (31) instead of Eq. (29) as a 
basic quantity is that the former (unlike the latter) is amenable 
to Pade summation in the test mass limit. 

5 As the energy is only a function of even powers of vq , the choice 
of using diagonal or subdiagonal (supcrdiagonal) is based on the 
order of v^ that is retained. For notational consistency, the 
indices on all approximants will refer to the power of vq . Other 
references define the indices on the energy approximants with 
respect to u^. 



approximants under variations of an (at the time) un- 
known coefficient in the 3PN center-of-mass energy, 
Damour, Jaranowski and Schafcr [24] found it more ro- 
bust to use the superdiagonal P-approximant instead of 
the subdiagonal P-approximant at 3PN order. 6 This sug- 
gestion was also adopted in Ref. [32] and will be used 
here; that is, we use subdiagonal P-approximants for 
1PN, diagonal for 2PN, and superdiagonal for 3PN. 

The P-approximants for the center-of-mass energy are 
defined as 
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\Ji + 4M - 1 
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where at 2PN order [21] 



M, 
(33) 



el(v n ) = -vh , 3 , t V ,J 35 \ { n , (34) 
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and at 3PN order [24] 
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2. Gravitational wave energy flux 

As originally pointed out in Refs. [47, 48], the flux 
function in the test-mass limit has a simple pole at the 
light-ring position (i.e., the last unstable circular orbit 
of a photon). Motivated by this, DIS introduced a new 
flux-type function 



fk(vci) = 1 - 



f'n 



Upole(") 



Fk{va;u) 



(37) 



with the suggestion that v po \ e be chosen to be at the light 
ring (pole singularity) of the new energy function. 

In order to construct well behaved approximants, DIS 
proposed to normalize the velocity vq entering the loga- 
rithms in Eq. (26) to some relevant scale which they chose 
to be wiso(' / ), where the last stable orbit (LSO) is defined 



Subdiagonal P-approximants were extended to 3PN order in 
Ref. [45], and LAL [46] software uses those P-approximants for 
the energy function. 



as the minimum of the energy. Also, they factored out 
the logarithms yielding 



h(vn) = f v 2 vh° 



x 1 
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>'n 
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Vpole(^) 
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(38) 



where £i and J-° & ~ ac are functions of Ti. Through 3.5PN 
order, £ 6 = -1712/105, £ 7 = 0, and jrj°s- fac = jr. with 
the replacement of vq — > v\ so in T§ [see Eq. (26)]. 

Finally, DIS proposed to define the P-approximant of 
the GW energy flux as 



F™M 
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1 - f'n/v po ie(^) 



fn(vn) 



(39) 



where 
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where P™[x] denotes Pade summation of the series x. 
DIS proposed to use the diagonal or subdiagonal P- 
approximants, depending on whether k = n + m is even 
or odd. Furthermore, DIS proposed to use v\ so (u) and 
w P oic(i / ) as the minimum and pole of the center-of-mass 
energy P-approximant of the same PN order. At 2PN 
(the order to which the PN expansion was known by 
DIS) u P oic is determined from the pole of the Pade energy 
function e|, yielding 



„2PN/ 



, 1 

L'polc \V) - -7= 



1 



1 



:±2, 
.if) 



(41) 



When the PN expansion was extended to 3PN order, it 
was found that none of the 3PN P-approximants have a 
physical pole. Therefore, somewhat arbitrarily, we will 
follow previous analyses and use the value (41) also at 
3PN order. We denote the P-approximants defined by 
Eqs. (39) and (33) as F^/E*. 

The denominator in the Pade summation of the GW 
energy flux can have zeros. They are called extraneous 
poles of the P-approximant [42] . It is desirable that these 
poles be located at high frequency (i.e., beyond the tran- 
sition from inspiral to plunge) . We shall see that depend- 
ing on the PN order and also the mass ratio, extraneous 
poles can be present at low frequencies. This could indi- 
cate poor convergence of the Pade summation. 

In Sees. IV B, VI B and VI C we shall investigate how to 
improve the closeness of the PN-approximants to the nu- 
merical data by varying 0,5 [16, 18, 26], v po \ c [18, 26] and 
also by introducing higher-order PN coefficients in the 



flux function. When varying i> po ic in the P-approximant 
at 3.5PN order, extraneous poles appear at low values of 
vq. Therefore, in order to push these poles to very high 
frequency, we follow the suggestion of Ref. [18], and use 
P-approximants at 4PN order, where the 4PN coefficient 
is set to its known value in the test-mass limit. This cure 
may fail for different mass ratios if new extraneous poles 
appear at low frequency. Furthermore the logarithm in 
the flux is not factored out as in Eq. (38), but treated as 
a constant when Pade summation is done. In this case 
the flux function is denoted F n . 

We notice that DIS motivated the introduction of the 
P-approximants first in the test-mass limit case by ob- 
serving much faster and monotonic convergence of the 
Pade energy, flux and waveforms with respect to Tay- 
lor energy, flux and waveforms. Quantitative tests of 
the convergence were done only for the Pade waveforms 
(see e.g., Tables III and IV in Ref. [21]), while for the 
flux and the energy conclusions were drawn qualitatively 
from Figs. 3 and 4 of Ref. [21]. DIS then conjectured 
that the comparable mass case is a smooth deformation 
of the test-mass limit case, and proposed to use close- 
to-diagonal P-approximants for the flux and the energy 
when v ^ 0. In the Appendix we perform a few conver- 
gence tests of the P-approximants of the flux function in 
the test-mass limit case, and conclude that whereas the 
P-approximants provide a better fit to the numerical flux 
at 5.5PN order, they do not accelerate the convergence 
of the Taylor series expansion of the energy flux. 



The EOB real Hamiltonian is 



iJ roal = M i / 1 + 2v 



H cS 



/' 



-M, (44) 



and we define H rcal = H lcai / \x. The T-approximants 
to the coefficients A(r) and D(r) in Eqs. (42) and (43) 
read [22, 24] 



fe+i 
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where 



Oo = l, ai = 2, a 2 = 0, a 3 (v)=2u, 

/ ^ f M 41 A 

a4(z/) = [t ~ 32* ) u > 

do = 1 , d\ = , diiy) = 6v, 
d 3 {v) = 2 (3v -26) v. 



(47) 



(48) 



In Sec. VI C, we will explore the flexibility of the EOB 
model by tuning the pseudo 4PN order coefficients a§(v) 
which we will take to have the following functional form 7 



a b (v) =a 5 v . 



(49) 



C. Non-adiabatic effective-one-body approximants 

The EOB model goes beyond the adiabatic approxima- 
tion and can incorporate deviations from the Keplerian 
law when the radial separation become smaller than the 
last stable circular orbit. 

Here we briefly review the main equations defining the 
EOB dynamics and refer the reader to previous papers 
for more details [15, 16, 18, 19, 22-24, 33]. The non- 
spinning EOB effective Hamiltonian is [22, 24]: 



i? cff (r,p) = ^H eS '(r,p) 



ix \A{r) 

1 



A(r) 
D(r) 



1 (n-p) 2 



2(4 - Zv) v (n ■ p) 4 



1/2 



(42) 



with r and p being the reduced dimensionlcss variables; 
n = r/V where we set r = |r|. In absence of spins the 
motion is constrained to a plane. Introducing polar co- 
ordinates (r, &,p r ,p<f,), the EOB effective metric reads 



In order to assure the presence of an horizon in the 
effective metric, we need to factor out a zero of A(r). 
This is obtained by applying the Pade summation [24]. 
Thus, the coefficients Ak(r) and D^(r) are replaced by 
the Pade approximants [24] 



A\{r) 



r (-4 + 2r + v) 
2r 2 + 2v + r v 



at 2PN order, and 



Mr) = 



Num(^) 
Den(4) ' 



(50) 



(51) 



with 



Num(A^) = r 2 [{a A {v) + 8v - 16) + r (8 - 2v)\ , (52) 

and 

Den(^) = r 3 (8 - 2v) + r 2 [ai(v) + Av] 

+r[2a 4 {v) + 8is]+A[is 2 + a i (iy)} 1 (53) 



dsi s = g^dx^dx v 



-A(r) c 2 dt 2 
-r 2 (d6> 9 



D(r) 



dr* 



Mr) 

sin 2 9 d<f) 2 ). (43) 



7 Note that what we denote as in this paper was denoted A in 
Ref. [16]. 
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at 3PN order. When exploring the flexibility of the EOB 
model, we use the following Pade approximant at 4 PN 
order [16, 26]: 



A\{r) 



Num(Aj) 
Den(AJ) ' 



(54) 



tangential velocity V$ = $r is set to V$ = wn = (& 1 / 3 , 
having assumed the Keplerian relation <£> 2 r 3 = 1 . It was 
then pointed out in Ref. [34] that the Keplerian relation 
becomes less and less accurate once the binary passes 
through the last stable orbit. A more appropriate ap- 
proximant to the flux would be 



with 



Num(A^) = r 3 [32-24^-4a 4 (*/)-a 5 (i/)] 

+r 4 [a 4 (f)-16 + 8f], (55) 



and 



Den(A\) = -a\{y) - 8a 5 (v) - % ai (v)v + 2o 5 (f)f - Vol? 
+r [-804(f) - 4o 5 (f) - 2a 4 (v)v - 16^ 2 ] 
+r 2 [-4ek(f) - 2a 5 {v) - 16f] 
+r 3 [-2a4(f)-a 5 (f)-8f] 
+r 4 [-16 + 04(f) +8f]. (56) 

For the coefficient D(r), the P-approximant used at 2PN, 
3PN, and 4PN order respectively are [16, 24, 26]: 



J$Cr)-l-£ 
D° 3 (r) 



r 3 + Qvr + 2^(26 - 3f ) 



(57) 
(58) 



Dl(r) 



r 4 + 6ur 2 + 2^(26 - 3f)r - d 4 (f) + 36f 2 ' 

(59) 



and we choose somewhat arbitrarily di{v) = 36f , so 
that D4 = D®. (We note that the value of d 4 does not 
affect much the EOB evolution [16].) The EOB Hamilton 
equations written in terms of the reduced quantities jj real 
and t = t/M, = fiM [23], are 



dr dH 



■ca.l 



■{r,Pr,PS>): 



dt dp r 

d$ ~ dH ieai . 

-~ = "=-5 {r,Pr,Ps>), 
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-(r,p r ,p 9 ). 



dt 



(60) 

(61) 

(62) 
(63) 



where for the <I> component of the radiation-reaction force 
a few approximants are available. Originally, Ref. [23] 
suggested the following Keplerian P-approximants to the 
flux 
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F n{ v n;v,v vo \ e ), 



(64) 



where F™ is given by the Pade flux in Eqs. (39) and 
(40). Here by Keplerian we mean that in the flux the 
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F™(V^;f,u po ie). 



(65) 



where V$ = <& rn . Notice that because the EOB Hamil- 
tonian is a deformation of the Schwarzschild Hamilto- 
nian, the exact Keplerian relation is $ 2 r^ = 1 with 
m = r [ip(r, p*)] 1 ' 3 and ip is defined following the ar- 
gument presented around Eq. (19) to (22) in Ref. [34]: 



1 

■0r 3 



= W?:, 



f dH(r,p r =0^^) 
V d P<t> 



(66) 

where w(r, p$) = A(r) ( 1 + -^ J . The value of p^ of cir- 
cular orbits are obtained by minimizing with respect to 
r the circular orbit Hamiltonian H(r, p r = 0,^) and it 
yields the following relation between r and p$ 



2j^4(r) 



r A I dr 



(67) 



By inserting Eq. (67) in the definition of i/>, and replacing 
all p^ except the one which implicitly appears in w(r^p^) 
we obtain 



i> 



l + 2??(y/w(r,p0)- 
r 2 dA(r)/dr/2 



(68) 



Finally, Refs. [32, 33] introduced another possible 
variation of the EOB flux approximants which use T- 
approximants for the flux given by Eq. (19), in either the 
Keplerian or non-Keplerian form, i.e. 



^T — 
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and 



11 K 



J n 
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F n {vci) , 



F n (V*). 



(69) 



(70) 



Note that the flux for the non-Keplerian EOB models 
are not simply functions of the orbital frequency Q. We 
denote the original E-approximants [22-24] which use the 
Pade flux (40) as F™/H p where H p is H lcal computed 
from A p and D p . Other E-approximants used in this 
paper are summarized in Table I. The initial conditions 
for Eqs. (60)-(63) are obtained following Ref. [23] and 
starting the evolution far apart to reduce the eccentricity 
to negligible values. 
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D. Waveforms 



The PN waveforms are obtained by substituting the 
orbital phase and frequency into the spherical harmonic 
mode (2,2) with amplitude corrections through 3PN or- 
der [49, 50] 
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For the adiabatic models, the orbital phase is obtained 
by rewriting the energy balance equation (12) as 



~dt 



F 



dE/dVL 



(72) 



and integrating this equation along with d^/dt = SI. 
The Taylor approximants are formed first by substituting 
F = F n and E = E n into Eq. (72). The P-approximant 
waveform is formed similarly by substituting F = F™ 
and E = E™ into Eq. (72). The TaylorTl and Pade 
approximants then numerically integrate Eq. (72). The 
TaylorT4 approximant is formed by first re-expanding 
the right side of Eq. (72) as a single Taylor expansion 
truncated at the appropriate order, and then numeri- 
cally integrating the resulting equation. The TaylorT2 
and TaylorT3 approximants perform the integration an- 
alytically. The various Taylor approximants are reviewed 
in Sec. HIE of Ref. [12]. 

For the non-adiabatic EOB models, the orbital phase 
is determined by solving Hamilton's equations (60)-(63). 

After computing /122, the appropriate time derivatives 
are taken to form /122 and ^| 2 . 



IV. COMPARISON WITH POST-NEWTONIAN 
APPROXIMANTS: ENERGY FLUX 

We now compare the numerical GW energy flux with 
predictions from PN theory. In Sec. IV A we present 
comparisons with T-, P- and E-approximants, and in 
Sec. IV B we explore ways of fitting the numerical flux 
by introducing higher-order PN coefficients and varying 
the value of v po i c away from v^^, [Eq. (41)]. 

The PN flux is derived as a function of frequency, so 
it is natural to perform this comparison as a function 
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FIG. 5: Ratio of GW frequencies uj and vj to orbital frequency, 
as a function of (twice) the orbital frequency, for different PN 
models. The GW frequencies lo and vo are defined in Eqs. (4) 
and (5). Solid lines correspond to 3.5PN, dashed and dotted 
lines to 3PN and 2.5PN, respectively. 



of frequency. One alternative, comparison as a function 
of time, would require computation of the PN phase as 
a function of time. This depends on the PN energy, so 
that a comparison with respect to time would mix effects 
due to flux and energy. Furthermore, comparisons with 
respect to time are sensitive to (and likely dominated by) 
secularly accumulating phase differences [51]. 

The PN flux is given in terms of the orbital frequency 
SI — see Eqs. (19) and (13) — so at first glance, it might 
seem natural to compare PN and NR energy fluxes at 
particular values of SI. However, the orbital frequency is 
gauge-dependent, and there is no simple relation between 
the NR orbital frequency and the PN orbital frequency. 
Nor is there a simple relation between the NR orbital fre- 
quency and any quantity measured at infinity (where the 
energy flux is defined). In particular, it is very difficult 
to determine the NR orbital frequency as a function of 
retarded time. In contrast, the frequency w (see Eq. (5)) 
of the GWs at infinity is an observable quantity, and is 
easily obtained from both PN formulae and from the NR 
simulation. Therefore, to achieve a meaningful compari- 
son, we compare the PN and NR energy flux at particular 
values of zu. 

In order to compute the PN flux as a function of zu, 
we need to find the mapping zupn : SI —> zu. In order 
to find this mapping, we must build a PN waveform as a 
function of SI and compute w as defined by Eq. (5). We 
construct the waveforms as described in Sec. HID. For 
the T-approximant of the flux, we will use the TaylorT4 
waveform. In Fig. 5 we plot both GW frequencies (de- 
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FIG. 6: Effect of choice of frequency. Shown are the PN fluxes 
for two representative PN-approximants, plotted (correctly) 
as function of vj and (incorrectly) as function of 20,. Plotting 
as a function of 20. changes the PN fluxes significantly relative 
to the numerical flux Fnr. 



fined in Eqs. (4) and (5)). We then invert the mapping to 
obtain J7pn = Wp^ : w — ► fl. So, given the PN flux -F(fi) 
from Sec. Ill, the flux as a function of the GW frequency 
is given by F(w) = F (ilpN (^)) ■ The relation f^p^n?) 
depends on the instantaneous evolution of the PN model 
around frequency 12, and is therefore (unfortunately) de- 
pendent on the PN model, in particular the choice of PN 
energy. This dependence, however, is local and will not 
lead to secularly accumulating differences. 

Notice from Fig. 5 that the orbital frequency and the 
GW frequency differ by ~ l%-3% at large frequencies, 
depending on the PN model and the PN order, and the 
difference in w between different PN models is about 5%. 
Because the energy flux is roughly proportional to tu 10 ' 3 
(more precisely, dlogF/dlog(A/tn) increases to ~ 3.6 at 
Mw = 0.15), the difference in the flux caused by using 
GW frequency from different PN models is about three 
to four times the difference in GW frequencies. Fig. 6 
illustrates this effect by intentionally plotting the PN flux 
versus the incorrect frequency Q. Because changing the 
PN model has a significant effect on the flux, we consider 
flux comparisons for several different PN models below. 

Note that for the flux comparison (and the compar- 
isons of the derivative of the energy in Sec. V), the PN 
waveforms are used only to define the mapping between 
f2 and w. The PN flux is taken directly from the PN flux 
expressions, e.g., Eq. (19), and not computed by apply- 
ing Eq. (3) to PN waveforms h(t). Equation (3) is used 
only to compute the numerical flux. 



A. Flux comparison 

Figure 7 plots the NR flux and the fluxes for the T-, 
P-, and E-approximants at 3.5PN order as a function 
of the GW frequency w computed from hi2- The T- 
approximant is TaylorT4 [12]. Along the top of this figure 
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FIG. 7: Comparison between the numerical energy flux and 
several PN approximants at 3.5PN order versus GW fre- 
quency w extracted from /122 in the equal-mass case. 
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FIG. 8: Comparison between the numerical energy flux and 
several PN approximants versus GW frequency zn extracted 
from /122 in the equal-mass case. We show the relative dif- 
ference between numerical flux and PN flux, as well as the 
estimated error of the numerical flux (blue bars, see Fig. 2). 
Solid lines represent 3.5PN models and NR; dashed and dot- 
ted lines correspond to 3PN and 2.5PN models, respectively. 
For notation see Table I and caption therein. 
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FIG. 9: Comparison of normalized energy flux F/Fj-j av ,t [see 
Eq. (73)] for the equal-mass case. Solid lines represent 3.5PN 
models and NR; dashed and dotted lines correspond to 3PN 
and 2.5PN models, respectively. For notation see Table I and 
caption therein. 



(as in several figures below) we indicate the number of 
gravitational wave cycles up to merger, where we define 
"merger" as the maximum of | "Fip | . Figure 8 zooms over 
the first 15 GW cycles. We notice that during the first 15 
GW cycles the numerical data are fit best by the P- and 
E-approximants at 3PN and 3.5PN order. At these low 
frequencies the NR flux is best matched by the Keplcrian 
and non-Keplerian EOB models and the Pade model. 

To more clearly show the behavior of the PN approx- 
imants, we plot in Fig. 9 the energy flux normalized by 
the Newtonian flux. The normalized flux is computed as 
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where for the same reason mentioned above, the New- 
tonian flux is expressed in terms of the GW frequency. 
Notice that the P-approximants and some of the E- 
approximants use the same Pade flux, but they start 
differing at Mm ^0.12 due to their different GW fre- 
quencies (obtained from an adiabatic and non-adiabatic 
evolution, respectively). The E-approximants with Ke- 
plerian and non-Keplerian flux increase less abruptly at 
high frequency than the P- and T-approximants. This 
is a consequence of non-adiabatic effects captured by the 
EOB model. Quite remarkably, the E-approximants with 
non-Keplerian fluxes are rather close to the NR result 
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FIG. 10: Cauchy convergence test of F/FNcwt for T- and 
P-approximants. We plot AF n+m = F n+m+1 — F n+m , and 
AF™ = F™ +1 - F™ for different values of vn- The T- and 
P-approximants are given by Eqs. (19) and (39), respectively. 
Note that the P-approximant has an extraneous pole at 1PN 
order at vn — 0.326. We use v iso = vf™, and v po i c — up 1 ™. 



for the entire range of frequency spanned by the sim- 
ulation. 8 We observe that somewhat accidentally the 
PN-approximants at 2.5PN order are also close to the 
numerical flux. 

The normalized NR flux starts to decrease at Mm ~ 
0.13. We notice that this behavior is rather different 
from the behavior of the normalized flux in the test-mass 
limit (see Fig. 19 in the Appendix). The E-approximants 
with non-Keplerian Pade or Taylor flux show a similar 
decreasing behavior at high frequency. 

Both Figs. 8 and 9 show that in the equal-mass case 
P-approximants fit the numerical results better than T- 
approximants. In numerical analysis, however, Pade 
summation is often used as a technique to accelerate the 
convergence of a slowly-converging Taylor series (e.g., see 
Tables 8.9 and 8.12 in Ref. [42]); hence it is natural to 
ask in the PN case whether Pade summation indeed ac- 
celerates the convergence of the series. In Table II we 
list the T- and P-approximants of F/Fj^ ewt computed at 
subsequent PN orders and for several values of vq [from 
left to right vq =0.1, 0.25 (i.e., beginning of the numeri- 
cal simulation), 0.3, 0.35, and 0.4.] In Fig. 10 we perform 



We notice that whereas the Keplcrian Pade-based (or Taylor- 
based) approximants to the flux differ from each other only when 
expressed in terms of the GW frequency, the non-Keplerian Pade- 
based (or Taylor-based) approximants to the flux differs from the 
others because their functional dependence on the frequency is 
different (e.g., compare Eq. (65) with Eq. (64)). 
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0.939 


-7.8434 
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1.091 


8.443 


1.5 


0.9681616 0.9686094 


0.9188 


0.9074 


0.940 


0.9069 


0.995 


0.924 


1.094 


0.967 


2.0 


0.9681512 0.9676191 


0.9184 


0.8850 


0.939 


0.8671 


0.993 


0.860 


1.091 


0.867 


2.5 


0.9675775 0.9676981 


0.8624 


0.8890 


0.799 


0.8754 


0.692 


0.875 


0.504 


0.893 


3.0 


0.9677265 0.9677247 


0.8951 


0.8914 


0.895 


0.8804 


0.928 


0.883 


1.022 


0.903 


3.5 


0.9677274 0.9677233 


0.8957 


0.8912 


0.897 


0.8798 


0.934 
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1.036 
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TABLE II: Normalized energy flux F/Fkcwt for the T- and P-approximants at subsequent PN orders for select velocities vn- 
vq — 0.25 corresponds to the start of the numerical simulation. The P-approximant flux is given by Eq. (39). Note that the 
P-approximant has an extraneous pole at 1PN order at va = 0.326. We use v\ BO = vf™ — 0.4456 and w po i c = ^p™ — 0.6907. 
We use boldface to indicate the range of significant figures that do not change with increasing PN order. 



a Cauchy convergence test and compute the difference 
between T- and P-approximants at subsequent PN or- 
ders. The figures do not suggest an acceleration of the 
convergence. We notice that in the equal-mass case P- 
approximants are converging more systematically than 
T-approximants. However, this fact seems to depend on 
the mass ratio, as can be seen by comparing Fig. 10 with 
Table IV and Fig. 21 in the Appendix which are obtained 
in the test-mass limit. 



On the fitting of the numerical relativity energy 
flux 



In view of building accurate analytical templates that 
can interpolate the NR waveforms during inspiral, merger 
and ringdown, we explore here the possibility of improv- 
ing the PN-approximants to the energy flux by introduc- 
ing phenomenological higher-order PN coefficients and/or 
by varying the value of w po ie ■ This study should be con- 
sidered a first exploration of the problem, demonstrating 
only the flexibility of the PN models. None of the quan- 
tities derived here should be used as the basis for further 
work. 

We will minimize the difference between the PN flux 
and the numerical flux by varying particular coefficients 
in the PN model. Ideally, the PN and numerical fluxes 
should be expressed as functions of w before taking this 
difference, so that the fluxes arc compared in a physi- 
cally meaningful way. Unfortunately, the calculation of 
w for the PN models is time-consuming, because for each 
trial value of the phenomenological coefficient it is neces- 
sary to compute a full waveform to determine the map- 
ping between w and £1. So instead, in this section we 
simply compare PN and numerical fluxes as functions of 
f2, where we define the numerical orbital frequency as 
f2 = tzj/2. In Fig. 6, we can see that the error intro- 
duced by the discrepancy between Q and zu/2 will be 
significant. As we will show in Sec. VI B, the waveforms 
produced using these "tuned" flux functions will improve 
agreement with the numerical waveform at a significant 
level. Nevertheless, the values derived in this section may 
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FIG. 11: Fitting several PN-approximants to the numerical 
flux. The a:-axis denote the orbital frequency f2. Because the 
numerical flux is computed as function of the GW frequency, 
we use for the numerical flux f2 = w/2. The blue bars indi- 
cate estimated errors on the numerical flux, see Fig. 2. For 
notation see Table I and caption therein. 



not be optimal. Thus, we emphasize that the results of 
this section constitute merely an exercise demonstrating 
the feasibility of adjusting the PN parameters to optimize 
the agreement of the PN flux function with numerical 
data. 

The least-squares fits are done on F (w) / F-^ cvit {w) [see 
Eq. (73)]. In the case of T-approximants, we fit for the 
unknown 4PN-order coefficient in Eq. (19) for the equal- 
mass case. We perform a least-squares fit of the 4PN- 
order function T$,(v = 1/4) = Ag + _Bgloguo over the 
orbital-frequency range MVt = 0.02-0.08 which starts af- 
ter the first 9 GW cycles. We obtain A 8 = -141 , B 8 = 
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102. We notice that when we perform the fit over the 
first 15 (or 20) GW cycles, spanning the frequency re- 
gion MCI = 0.0168-0.0235 (Mfl = 0.0168-0.0283), the 
agreement becomes worse. The resulting flux is shown in 
Fig. 11. The relative difference with the numerical flux 
is at most ~ 0.8%. 

We repeat this analysis in the case of P-approximants. 
Because the latter also depend upon i> po ie, we perform 
two least-squares fits. In the first fit, we fix v po i c to the 
value given by Eq. (41) and apply the least-squares fit to 
fg(v = 1/4) obtaining A 8 = -1382 , B 8 = 197. 

In the second fit, we vary v po i . When varying v po \ e 
in the P-approximant at 3.5PN order, extraneous poles 
appear at low values of vn- Therefore, in order to push 
these poles to very high frequency, we follow the sugges- 
tion of Ref. [18], and use P-approximants at 4PN order, 
where the 4PN coefficient is set to its known value in the 
test-mass limit. Furthermore the logarithm in the flux is 
not factored out, but treated as a constant when Padc 
summation is done. This cure may fail for different mass 
ratios if new extraneous poles appear at low frequency. 
The least-squares fit gives v po i = 0.74. All the results 
for the P-approximants are displayed in Fig. 11, where 
we also show the T- and P-approximants at 3.5PN order 
without any fit. 

Figure 11 might suggest that by introducing higher- 
order PN coefficients in the flux, the numerical flux can 
be fit better by T-approximants than by P-approximants. 
However, this result can depend on the use of orbital 
frequency instead of GW frequency. In Sec. VIC (see 
Fig. 18) we employ the fit values obtained in this study 
and show phase differences between NR and tuned EOB 
models. 

Finally, we attempted to extract PN coefficients higher 
than 3.5PN order from the numerical flux, as was done at 
2PN, 2.5PN and 3PN order in Ref. [48] in the test-mass 
limit. Unfortunately, the differences between numerical 
flux and T-approximants are so large — even at the begin- 
ning of the numerical waveform — that we were not able 
to extract even known PN coefficients, like the ones at 
3PN and 3.5PN order. Thus, to fit unknown PN coef- 
ficients would require a numerical simulation with more 
cycles starting at lower frequency. 



V. ESTIMATION OF THE (DERIVATIVE OF 
THE) CENTER-OF-MASS ENERGY 
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FIG. 12: GW frequency derivative vj for the numerical rela- 
tivity simulation and various PN approximants at 3.5PN or- 
der. For notation see Table I and caption therein. 



energy, as can be seen by rewriting the energy balance, 
Eq. (12), in the form 



dvj 
~~dt 



F 



dE/dvo 



(74) 



Therefore, we begin this section with a comparison be- 
tween numerical w and the predictions of various PN- 
approximants. For the PN-approximants, we compute 
/122 as usual (i.e., using energy balance to compute the or- 
bital frequency derivative fl), and take a time derivative 
to obtain hii and extract vj from it. The waveform /122 
for the E-approximants is computed using Eqs. (42), (44), 
(45) and (46) in Sec. IIIC. Figure 12 plots the numerical 
vj and its value for T-, P- and also E-approximants at 
3.5PN order. 

In order to emphasize differences between the different 
vj, we normalize the data in Fig. 12 by the Newtonian 
value of vj, 



In the previous section, wc analyzed and compared PN 
and numerical energy fluxes. The energy of the binary is 
the second fundamental ingredient in the construction of 
adiabatic PN-approximants. Unfortunately, there is no 
way to extract the energy for the numerical simulation 
as a function of a gauge-invariant quantity such as the 
GW frequency, so that it is impossible to compare PN 
and NR energies directly. The frequency derivative, 137, 
however, is easily accessible in the numerical data, and, in 
the adiabatic approximation is intimately related to the 
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The normalization is used only to eliminate the leading- 
order behavior of the various curves in Fig. 12; therefore, 
to compute the denominator of Eq. (75) we have sim- 
ply substituted zu/2 rather than Q into the Newtonian 
formula for the frequency derivative. 

The normalized frequency derivatives are shown in 
Fig. 13. At low frequencies, vj is very challenging to 
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FIG. 13: Comparison of zu for the numerical results and var- 
ious PN approximants. Dotted, dashed and solid lines corre- 
spond to 1.5PN, 2.5PN and 3.5PN models, respectively. For 
notation see Table I and caption therein. 



FIG. 14: Comparison of PN w with a heavily smoothed ver- 
sion of the numerical vu. Solid lines represent 3.5PN mod- 
els and NR; dashed and dotted lines correspond to 3PN and 
2.5PN models, respectively. For notation see Table I and cap- 
tion therein. 



compute in numerical simulations, resulting in compar- 
atively large numerical uncertainties. Therefore, for fre- 
quencies Mvu < 0.045 we can merely conclude that PN 
and NR are consistent with each other (i.e., are within 
the numerical error bars of about 10 per cent). 

The 3.5PN Taylor T4 model (labeled F 7 /E 6 T4) agrees 
very well with the numerical simulation up to Mvu w 0.1; 
this observation is consistent with the excellent agree- 
ment between TaylorT4 (3.5PN) and the numerical sim- 
ulation observed in Boyle et al. [12], who compared up to 
this frequency. Beyond mvu = 0.1, however, vu/vu-^ ewt for 
Taylor T4 continues to increase (as for all other Taylor 
and Pade models considered here), whereas for the nu- 
merical simulation, vu/vu-Newt flattens (this behavior was 
also observed in Ref. [18].) Only the E-approximants 
at 3.5PN order reproduce the flattening of vu/vu-Ncwt at 
high frequencies, with the closest being the one which 
uses the non-Kcplcrian Pade flux ( nK F|). Because the 
frequency derivative is the relevant quantity that deter- 
mines the phase evolution, the turning over of vu/vuNewt 
for the non-adiabatic models in Fig. 13 suggests that, at 
high frequency, non-adiabatic analytical models might be 
superior to adiabatic models. 

If sufficient smoothing is applied to the numerical vu it 
becomes a smooth curve even at low frequencies. Fig- 
ure 14 presents a comparison between such a heavily 
smoothed numerical curve and the PN-approximants. As 
already pointed out, all PN approximants are consistent 
to within our estimated numerical errors at low frequen- 
cies. However, the NR result in Fig. 14 is notably closer 



to the 3.5PN approximants than to lower order PN ap- 
proximants. This good agreement provides a further val- 
idation of the numerical code used in Boyle et al. [12]. It 
also indicates that our error analysis in Sec. II may be 
overly conservative. 

Our comparisons of vu reveal a lot of information about 
the PN approximants. However, vu depends on both flux 
and energy (see Eq. (74)), and so these comparisons do 
not yield information about flux or energy separately. To 
isolate effects due to the PN energy, we rearrange Eq. (74) 
further, such that it yields in the adiabatic approxima- 
tion the derivative of the center-of-mass energy for the 
numerical simulation: 



dE 

dvu 
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NR 
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NR 



(76) 



The relative error in [dE /dvu\^ K is obtained as the root- 
square-sum of the relative errors of flux and frequency 
derivative (see Figs. 2 and 4). In Fig. 15 wc compare 
the latter with T-, P- and E-approximants. For adia- 
batic T4 and Pade models, we compute dE/dvu by tak- 
ing derivatives of E(il) in Eq. (14) with respect to fi 
and then expressing the derivative in terms of vu(£l). For 
non-adiabatic EOB models, we compute dE/dvu from the 
ratio of Fpn and [dvu/dt]p^ as obtained from Figs. 7 and 
12. The closeness between the numerical result and adi- 
abatic PN-approximants is expected only in the range 
of frequencies over which the balance equation and the 
adiabatic approximation are valid. The upper panel of 
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FIG. 15: We compare dE/dzu versus GW frequency w for nu- 
merical relativity [see Eq. (76)] and PN approximants. Solid 
lines represent 3.5PN models and NR; dashed and dotted lines 
correspond to 3PN and 2.5PN models, respectively. For no- 
tation see Table I and caption therein. 



Fig. 15 shows the Taylor and Pade adiabatic models. The 
plot suggests that around Mm ~ 0.08 non-adiabatic ef- 
fects are no longer negligible. At lower frequencies, both 
3.5PN order adiabatic approximants (Pade and Taylor 
T4) match the numerical result very well. Taylor T4 at 
2.5PN matches well, too, although its frequency deriva- 
tive w and flux differ significantly from NR (see Figs. 13 
and 9). The T-approximant at 3.5PN order is closest to 
the numerical result. The lower panel of Fig. 15 shows 
the non-adiabatic E-approximants. We notice that the 
non-adiabatic models, especially at 3.5PN order, follow 
quite nicely the behavior of the numerical derivative of 
the center of mass energy. The E-approximant with non- 
Keplcrian flux is closest to the numerical result. This 
analysis emphasizes again the relevance of including non- 
adiabatic effects in the analytical model [23] . 



VI. COMPARING WAVEFORMS 

Here we compare the numerical waveform to various 
PN waveforms, basically extending the analysis of Boyle 
ct al. [12] to include Pade and EOB waveforms. Because 
the (2, 2) mode dominates the waveform for an equal- 
mass non-spinning binary, we restrict the comparison to 
only this mode. As in [12], we use 'J'! 2 and the GW 
phase and frequency u defined by Eq. (4) when compar- 
ing waveforms. 

For the comparisons presented in this section, the 



uncertainty in the phase of the numerical waveform is 
roughly 0.02 radians. This number includes numerical 
errors (e.g. due to convergence and extrapolation of the 
waveform to infinite extraction radius), as well as mod- 
elling errors due to slightly nonzero eccentricity and spin 
of the numerical simulation; see Ref. [12] Sec. V. for de- 
tails. We note that the modelling errors have decreased 
since the analysis in Ref. [12] because the new match- 
ing procedure reduces the impact of eccentricity, and be- 
cause the more sophisticated spin-diagnostics presented 
in Ref. [52]) resulted in a smaller bound on the residual 
spin. 



A. Matching procedure 

Each PN waveform has an arbitrary time offset, to, and 
phase offset, (fio with respect to the NR waveform. The 
procedure used by Boyle et al. [12] — as well as in vari- 
ous other papers before it, such as [10, 11] — sets these 
constants by ensuring that the GW phase and frequency 
match at a fiducial time. Unfortunately, when matching 
at low frequency this method is sensitive to noise and 
to residual eccentricity in the numerical waveform, and 
does not easily translate into a robust and automatic 
algorithm. Since we want to match as early as possi- 
ble (where we expect the PN approximants to be valid), 
we propose to use, instead, a matching procedure which 
achieves the same goal, but extends over a range of data. 
This procedure is similar to the one proposed by Ajith 
ct al. [17], but whereas we match only the GW phase, 
Ajith et al. match the entire gravitational waveform — 
including the amplitude — and include an overall ampli- 
tude scaling. This method can be easily implemented as 
a fairly automatic algorithm, robust against noise and 
residual eccentricity. 

Using the phase of the numerical and PN waveforms, 
we define the quantity 



S(At,A^>) 



[(fan® - <fa N (t - At) - Acj)] 2 dt . 



. (77) 
Here, t\ and £2 represent the chosen range over which to 
compare. Minimizing this quantity by varying the time 
and phase offsets At and A(f> produces the optimal val- 
ues for these quantities in a least-squares sense. Then 
to compare PN and NR waveforms, we compare the (un- 
changed) NR waveform with an offset PN waveform de- 
fined by 



<J? 



4,PN 



(t) =A PN (t + At)e 



-i[0P N (t+At)+A0] 



(78) 



With reasonable first guesses for At and A<p, the func- 
tion S is quite nicely paraboloidal. Thus, even sim- 
ple minimization routines work well. However, in cases 
where speed is an issue, the problem can be reduced to 
one dimension. For a given value of At, the optimization 
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over A<f> may be done analytically by setting 



A0(A£) = 



It" [<faa.(t) - <t>p N (t - At)] dt 



ti-U 



(79) 



Using this value of Acf> for a given value of At decreases 
the number of function evaluations needed to find the 
minimum. This can be very useful for large data sets, or 
situations where many such matches need to be done. 

The choice of £1 and £2 involves some degree of judg- 
ment. Preferably, t\ should be as early as possible, while 
not being contaminated by junk radiation. We choose 
ti = 1100M, corresponding to Ml) = 0.037. Similarly, 
£2 should be as early as possible, but far enough from 
£1 so that the integration averages over the noise. In 
addition, the effects of the small but nonzero orbital ec- 
centricity show up as oscillations in the phase, as can 
be seen, for example, in the range £ G [1100, 1900]M in 
Fig. 17. We would like £2 to be large enough so that 
the integration averages over several cycles of this oscil- 
lation, thus resulting in less bias due to eccentricity. Here 
we use £ 2 = 1900M, corresponding to Muj = 0.042. We 
have checked that changing the values of £1 and £2 by 
±100M changes the resulting phases by less than a few 
thousandths of a radian through the end of the numerical 
waveform. 

This method is quite similar to the one suggested in 
Ref. [17]. However, here we consider only the phase and 
not the amplitude of the waveform. Because we restrict 
the analysis only to the (2, 2) waveform mode of an equal- 
mass binary and compare only the phase and not the 
amplitude, we think it is reasonable to have neglected 
the amplitude in the matching procedure. 



B. Pade waveforms 

In Fig. 16 we plot the phase difference between the 
numerical, T- and P-approximants [21, 24, 32] at the 
times when the numerical waveform reaches GW fre- 
quencies Muj = 0.063 and Muj = 0.1. The phase dif- 
ferences are plotted versus the PN order. The phase dif- 
ference at Muj = 0.1 of the P-approximant at 3.5PN 
order is —0.12 radians. When comparing with generic 
Taylor approximants, we notice that the phase differ- 
ences of the P-approximants are less scattered as the 
PN order is increased. This might be due to the fact 
that P-approximants of the energy flux are closer to the 
NR flux, especially for lower vn where the phase accu- 
mulates the most. Figure 16 could be contrasted with 
Tables III and IV of Ref. [21] which show the overlaps 
between the numerical waveform and P-approximants at 
subsequent PN orders, in the test-mass limit case. The 
behavior of the P-approximants in Fig. 16 are consistent 
with the behavior of w seen in Fig. 13: At 1.5PN, Pade 
has smaller w than the numerical simulation, at 2.5PN, 
Pade has larger w. Consequently, $pn — $nr is negative 
at 1.5PN order and positive at 2.5PN order. For 3.5PN 



a 
a 

'-6 



Z 

0- 

-©- 



.5 



1 0.063 

* 


1 I Q Q 


t=t * 

0.1 

* * 
° 

* 

* * 


• 


i ' " " ° 

* H * 



9 

* 5 

WW ♦ 

1,1,1,1, 


.5 

-5 


♦ 
O 


* E-approximant 

* P-approximant 
O TaylorTl 

x TaylorT2 
□ TaylorT3 

* TaylorT4 




, 






1 2 3 
PN order 



1 2 3 
PN order 



FIG. 16: Phase differences between the numerical wave- 
form, and untuned, original EOB, untuned Pade, and Tay- 
lor waveforms, at two selected times close to merger. The 
E-approximants are F™/H p , while the P-approximants are 
F™/Ep (see Table I and caption therein). Waveforms are 
matched with the procedure described in Sec. VI A and phase 
differences are computed at the time when the numerical sim- 
ulation reaches Mui = 0.063 (left panel) and Mu> = 0.1 (right 
panel). Differences are plotted versus PN order. Note that at 
1PN order the Pade flux has an extraneous pole at v — 0.326 
causing a very large phase difference. The thick black line 
indicates the uncertainty of the comparison as discussed in 
Sec. VI, |$pn - $nr| < 0.02 radians. 



order, the P-approximant in Fig. 13 agrees very well with 
the numerical simulation (at least for Mm < 0.1), which 
translates into excellent agreement in Fig. 16. 

In Fig. 17 we explore the possibility of reducing the 
phase differences between the numerical waveform and 
P-approximants: By (i) varying v po \ e or (ii) introducing 
the pseudo 4PN order coefficient J-%{y = 1/4) = Ag + 
Bg, log vn in the energy flux. We tune the coefficients by 
minimizing the sum of the squares of the phase difference 
a t io.063 an d £0.1 ■ We find that if v po i c = 0.633, the 
P-approximant F^/E^ has a maximum phase difference 
before Muj = 0.1 smaller than the numerical error in the 
simulation. A similar result is obtained for the the P- 



approximant pF^/E^ if we use u po ic 



and tune Ax 



-493, B 8 = 330. 



£E 



0.6907, 



C. Effective-one-body waveforms 

In Fig. 16 we also plot the phase differences be- 
tween the numerical and the untuned, original E- 
approximants [22-24] F™/H p . At 3.5PN order the phase 
difference at Muj = 0.1 is 0.50 radians. We also computed 
the phase differences at Muj = 0.1 of the E-approximants 
nK Fi/H 7 , nK F 7 /H 7 and F 7 /H 7 and found 0.45, 2.56 and 
2.7 radians, respectively. Thus, for untuned EOB models 
it is crucial to have introduced the Pade flux. When con- 
trasting the original E-approximants with generic Taylor 
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FIG. 17: Phase differences between untuned and tuned P- 
approximants and NR. The untuned P-approximant is F'i jE\ 
(wiso = fi 2 S o N ! u poic = Wpole)- The tuned P-approximants are 



F4JE2 and tunable v po \ c (v\ B , 



Vi s 



, "pole = Wp i e ) with tunable At 



') and pF\/Ei (i; lso = 
and Bs- In all cases, 



waveforms are matched over t — r* £ [1100, 1900]M. 



approximants, we find that the phase differences are less 
scattered as the PN order is increased. However, de- 
spite the fact that the Pade-based EOB flux is closer to 
the numerical flux (see Figs. 8 and 9), untuned, original 
E-approximants accumulate more phase difference than 
P-approximants. This could be a consequence of the fact 
that independently of the flux and the energy functions, 
what seems to matter is the way the equations of motions 
are solved to get the phasing. 

Because of the reduction of the dynamics to a few cru- 
cial functions determining the inspiral evolution [22, 23, 
25], notably A, D and J 7 , and because of the rather sim- 
ple procedure to match the inspiral(-plunge) waveform to 
the ringdown waveform, the EOB model turned out to 
be particularly suitable for matching the full numerical 
waveforms [9, 16, 18, 20, 27]. In view of a future study 
which will include merger and ringdown, we start here ex- 
ploring the possibility of improving the agreement with 
numerical waveforms by tuning the pseudo 4PN order 
coefficients 05, As and B% and/or, if present, the pole lo- 
cation v po \ c . In the lower panel of Fig. 18, using different 
i'p°ic values, we show the phase differences computed at 
^0.063 an d to.i as functions of the unknown PN-expansion 
coefficient 05 [see Eq. (49)]. As first pointed out and dis- 
cussed in Ref. [18] (see e.g., Fig. 3 therein), we find that 
there is a strong degeneracy between 05 and u po ic- In 
fact, for different v po i e values, the curves in Fig. 18 are 
almost identical except for a shift in as . Although in this 
test we use the E-approximant F^ /pHg,(v\ so = i'f S o N ), we 
find that this degeneracy appears in all E-approximants 
considered. 
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FIG. 18: The upper panel shows phase differences versus 
time (lower a;-axis) and versus GW frequency (upper i-axis) 
for several tuned and untuned E-approximants. For the tuned 
models, the optimal 0,5 and v po \ e values displayed in Table III. 
In the lower panel we show phase differences between nu- 
merical and E-approximants computed at to. 063, £0.1, and the 
end of the numerical simulation £0.16, as functions of as. For 
the same color and style, the curve with the steepest slope 
corresponds to £0.16 and the curve with the smallest slope 
corresponds to £0.003 (For notation see Table I and caption 
therein). 



To obtain the optimal a§ and v po \ that minimize phase 
differences during the entire numerical simulation, we 
first choose an arbitrary u po i e in the range of degeneracy. 
Then, we determine the 0,5 value by minimizing the sum 
of the squares of the phase difference at to.063 and £0.1 • 
In the upper panel of Fig. 18, we show phase differences 
in time and GW frequency for several E-approximants 
using those optimal a§ and v po \ values, which arc given 
in Tabic III. In Fig. 18, we also show phase differences 
for E-approximants with pseudo 4PN order coefficients 
determined by the flux fit of Sec. IV B (see Fig. 11) and 
tunable 05. The optimal 0,5 values are shown in Tabic III. 
The smaller phase differences along the entire inspiral 
are obtained with the E-approximants with Pade flux 
F£/pH 8 (vuo = uf™) and tunable u po ic,a5 and Taylor 
flux pFg,/pHg, with tunable As, B$, 05. We notice that for 
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EOB model and fixed parameters 


as 


vpole 


nK Ft/pH 8 




29.78 


0.52 


Ft/ P H S 


„ 2PN 
Ulso — Vl BO 


39.35 


0.55 


pFs/pHg 


A 8 = -141, B 8 = 102 


5.32 


N/A 


pFi/pHs 


,4 8 = -1382,B 8 = 197, 

2PN ,, 2PN 


-3.10 


N/A 



TABLE III: Optimal as and w po i c that minimize phase differ- 
ences between tuned EOB models and the numerical simula- 
tion. 



t > io.i the phase difference increases more abruptly for 
the latter model. In the best case, the absolute phase dif- 
ference during the entire numerical simulation is within 
the numerical error, i.e., within 0.02 radians. The choice 
of the best tuned E-approximant [15, 16, 18-20] will be 
determined once merger and ringdown are included, and 
when long and accurate comparisons with numerical sim- 
ulations are extended to BBH with mass ratio different 
from one. 

Finally, in Ref. [18], Damour and Nagar extracted the 
data of the numerical simulation used in the present pa- 
per from one of the figures of Ref. [12] and compared 
those data with the EOB approach. They found for their 
"non-tuned" EOB model phase differences ±0.05 radi- 
ans. This phase difference is smaller than the phase dif- 
ferences we discuss in this paper for untuned EOB mod- 
els (see Fig. 16 and discussion around it). However, we 
notice that ±0.05 radians in Ref. [18] refers to half the 
maximum phase difference accumulated over the entire 
evolution when matching the numerical and EOB phases 
at Moo — 0.1. By contrast, in this paper, and in par- 
ticular in Fig. 16, we match numerical and EOB phases 
in a time interval and compute the phase differences at 
Ma; = 0.1. 

Moreover, we observe that their "non-tuned" EOB 
model is not really untuned, because it uses the Pade 
summation of the radial potential at 4PN order and then 
sets 05 = 0. This is not equivalent to using the radial 
potential at 3.5PN order with 05 = 0. In fact, to recover 
the 3.5PN order Pade radial potential from the 4PN order 
Pade potential one should use 05 = —17.16. They also 
use the non-Keplerian flux at 4PN order F 4 which is 
different from the 3.5PN order one nK F|. For our un- 
tuned EOB model at 3.5PN order which uses nK F% and 
the EOB dynamics at 3PN order, if wc apply Ref. [18] 
procedure and compute half the maximum phase differ- 
ence when matching the numerical and EOB phases at 
Mui = 0.1, wc find a phase difference of ±0.18 radians 



energy flux, and GW frequency derivative. Imposing 
the balance equation, we also estimate the (derivative 
of) ccnter-of-mass energy. We compare these quanti- 
ties to those computed using adiabatic Taylor T4 and 
Pade [21, 24, 32], and non-adiabatic EOB PN approx- 
imants [22-24]. 

We find that for the first 15 GW cycles, the 3.5PN 
order T-approximant and the 3.5PN order untuned P- 
and E-approximants (sec Table I) reproduce the numeri- 
cal results for energy flux, GW frequency derivative and 
(derivative of) center-of-mass energy quite well (see Figs. 
8, 9, 13, 14, and 15), but with interesting differences. 

We attempted to study the convergence of the PN ex- 
pansion for the energy flux. 9 We find that Pade approxi- 
mants to the flux introduced in Ref. [21] do not accelerate 
the convergence of the Taylor series, but are closer to the 
numerical flux than are the T-approximants. In particu- 
lar, the Taylor flux at all orders through 3.5 PN is outside 
the numerical flux error bars even ~ 25 GW cycles be- 
fore merger (see Fig. 8). We find that the non-adiabatic 
non-Keplerian E-approximants to the flux at 3.5PN or- 
der are within ~ 2% of the numerical flux over the entire 
frequency range we consider (see Fig. 9). 

Quite interestingly, in the equal-mass case the numer- 
ical normalized energy flux F/F^ cwt starts decreasing at 
high frequency during the late part of the inspiral and 
blurred plunge (see Fig. 9). This differs from the be- 
havior of F/F-Ncwt in the test-mass limit (see Fig. 19). 
Both the Taylor and Pade-based E-approximants with 
non-Keplerian flux [34] show a similar decreasing behav- 
ior at high frequency. This fact suggests that if a pole 
is present in the energy flux of equal-mass binaries, it is 
located at a larger frequency than that at which the com- 
mon apparent horizon forms. As seen in Sec. IV B, when 
fitting for Wp i c we obtain v po \ e {v = 1/4) = 0.74, which is 
to be contrasted with the test-mass case v po \ c {v = 0) = 
1/a/3 rs 0.58. These values of v vo \ a correspond to orbital 
frequencies MQ = 0.405 and Mfl = 0.192, respectively. 

For the GW frequency derivative w, we find that at low 
frequency the Taylor, Pade and EOB models at 3.5PN or- 
der are within the numerical error (see Fig. 13). At high 
frequency, as already observed in Ref. [18], only the non- 
adiabatic E-approximant has a GW frequency derivative 
that flattens out, as does the numerical result. The non- 
Keplerian E-approximant at 3.5PN order is closest to the 
numerical data (see Fig. 14). 

When estimating the derivative of center-of-mass en- 
ergy dE /dm, we expect the numerical result and adi- 
abatic PN-approximants to be close only in the range 
of frequencies over which the balance equation and the 
adiabatic approximation are valid. We find that this 



VII. CONCLUSIONS 

In this paper, using a highly accurate and long numer- 
ical simulation [12] of a non-spinning equal-mass black 
hole binary, we compute the gravitational waveform, GW 



9 We also tried to apply the criterion suggested in Ref. [53] to 
assess the region of validity of the PN series for the flux in the 
equal-mass case. Unfortunately, the numerical simulation starts 
at too high a frequency, when the Taylor series at 3.5PN order 
seems to already be outside the region of validity. 
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range of frequencies is Mw < 0.08 (see Fig. 15) for the 
2.5PN T-approximant and all the 3.5PN approximants. 10 
At higher frequency, the 3.5PN order non-adiabatic E- 
approximants are closer to the numerical dE/dw than 
are the adiabatic approximants, and the non-Keplerian 
E-approximant is the closest. 

Applying a new matching procedure, we compared the 
numerical waveforms with TaylorT4, Pade, and EOB 
waveforms. We find that the accumulated phase dif- 
ference from the numerical solution at Mu> = 0.1 is 
—0.12 radians for the untuned 3.5PN P-approximant [21, 
24, 32], 0.50 radians for the untuned, original 3.5PN E- 
approximant [22-24], and 0.45 radians for the untuned 
non-Keplerian [34] 3.5PN E-approximant (see Fig. 16). 
Although those phase differences are larger than for 
3.5PN TaylorT4 (-0.04 radians), the phase differences 
for the P-approximants are less scattered as a function 
of PN order than are the phase differences for generic 
Taylor approximants. 

The analyses of the flux, GW frequency derivative and 
(derivative of the) center-of-mass energy emphasize again 
the importance of including non-adiabatic effects dur- 
ing the last stages of inspiral [23] . Roughly, we can say 
that non-adiabatic effects are no longer negligible start- 
ing from a frequency Mw ~ 0.08-0.12, as can be seen 
in Figs. 9, 13, and 15. As seen in these figures, non- 
adiabatic E-approximants can capture some of the rele- 
vant features of the late time evolution. We expect that 
by further improving these models by fitting higher-order 
PN coefficients to the numerical data, they will become 
excellent candidates for developing an analytic template 
bank of coalescing BBHs [9, 16, 18, 20, 27]. 

In this paper we started to explore the possibility of 
reducing the phase differences between numerical and E- 
approximant waveforms by fitting the unknown parame- 
ters as, F%, and Wp ic (sec Fig. 18). As a first step, for 
several E-approximants we searched for a local minimal 
phase difference by varying 05, F%, and w po ie- We found 
that we were able to reduce phase differences to below 
the numerical uncertainty. In a future work which will 
include merger and ringdown, we plan to determine the 
region of the parameter space (05, T$, u po ic) in which 
the phase difference is within the numerical uncertainty 
of the simulation. 
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APPENDIX: PADE APPROXIMANTS TO THE 
ENERGY FLUX IN THE TEST PARTICLE LIMIT 

In the test-mass-limit case the GW energy flux is 
known through 5.5PN order [44]. The explicit coefficients 
entering Eq. (19) for i > 8 and v = can be read from 
Eqs. (4.1) and (4.2) of Ref. [21]. 

In Fig. 19 we compare the normalized energy flux func- 
tion [43] F'/Fkowt to the T- and P-approximants. To eas- 
ily compare Fig. 19 with the other figures in the paper, 
we plot quantities as functions of the approximate GW 
frequency defined by 2MfL As noticed in Ref. [21], the 
P-approximants approach the numerical data more sys- 
tematically The differences between different PN orders 
are difficult to see in Fig. 19. To obtain a clearer view, 
Fig. 20 plots the differences between PN flux and nu- 
merical flux at four fixed frequencies. Fig. 20 shows this 
somewhat better behavior of Pade; however, the Padc- 
approximants show little improvement between PN or- 
ders 3.5 and 4.5, and at order 5 there occurs an extra- 
neous pole. At frequency 2Mfl = 0.04, P-approximants 
with order > 2.5 are within 0.5 percent of the numeri- 
cal data, as are T-approximants with order > 3.5. Good 
agreement at low frequency is rather important because 
that is where the majority of the waveform phasing ac- 
cumulates. 

Table IV and Fig. 21 test the internal convergence of 
T- and P-approximants without referring to a numeri- 
cal result. Table IV displays the flux at all known PN- 
orders at select frequencies, with boldface highlighting 
the digits that have already converged. Although the 
Pade summation does not accelerate the convergence, the 
P-approximant at 5.5PN order is closest to the numerical 
data (see Fig. 20). 

Comparing Table IV with Table II, and Fig. 21 with 
Fig. 10 we observe that the P-approximants converge 
more systematically in the equal-mass case than in the 
test- mass limit. This is also evident by comparing Fig. 20 
with Fig. 8: We see that P-approximants at 3PN and 
3.5PN orders are inside the numerical flux error whereas 
T-approximants at all orders through 3.5 PN arc outside 
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PN order 


v n = 0.1; 2M Q = 0.002 


v n = 0.3; 2MQ, = 0.054 


vn = 0.4; 


2MQ = 0.128 


(n+m)/2 


-T 7i-\-m,/ ^Ncwt 


-FnV-pNcwt 


-T n-\-m / ^Ncwt 


F™ /Fviewt 


-F?i + m/-fNcwt F." l /FNcwt 


0.0 


1.0000000000 


1.20948977 


1.0000 


2.0817 


1.000 


3.255 


0.5 


1.0000000000 


1.03092783 


1.0000 


1.3699 


1.000 


1.923 


1.0 


0.9628869047 


0.94287089 


0.6660 


-0.9467 


0.406 


-12.52 


1.5 


0.9754532753 


0.97587569 


1.0053 


0.9916 


1.210 


1.201 


2.0 


0.9749604292 


0. 97462770 


0.9653 


0.9337 


1.084 


1.031 


2.5 


0.9745775009 


0.97469475 


0.8723 


0.9422 


0.692 


1.063 


3.0 


0.9747307757 


0.97471937 


0.9710 


0.9465 


1.227 


1.069 


3.5 


0.9747206248 


0.97471854 


0.9488 


0.9460 


1.061 


1.066 


4.0 


0.9747182352 


0.97471874 


0.9369 


0.9462 


0.952 


1.067 


4.5 


0.9747194262 


0.97471859 


0.9559 


0.9461 


1.190 


1.066 


5.0 


0.9747192776 


0. 97471930 


0.9479 


1.1178 


1.051 


1.037 


5.5 


0.9747192763 


0.97471928 


0.9485 


0.9493 


1.073 


1.091 



TABLE IV: Normalized energy flux F/FNewt in the test-mass limit for the T- and P-approximants at different PN orders and 
at three different frequencies. We use boldface to indicate the range of significant figures that do not change with increasing 
PN order. 



the numerical flux error bars even ^25 GW cycles be- 
fore merger. However, as the Padc approximant does not 
converge faster, it is not immediately clear whether sim- 



ilar superior behavior of Pade can be expected for more 
generic binary black holes. 
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